Part 2 Fire Selection
2.1 Set Up
2.1.1 Libraries
library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(gridExtra)
library(knitr)
library(rasterVis)
library(RColorBrewer)
library(spData)
library(forcats)
library(cowplot)
library(rgeos)2.1.2 USDA National Forest Type Group Dataset
Conifer Forest Type Groups: Douglas-Fir, Fir-Spruce-Mountain Hemlock, Lodgepole Pine
# forest type groups and key
conus_forestgroup <- raster('data/forest_type/conus_forestgroup.tif')
forest_codes <- read_csv('data/forest_type/forestgroupcodes.csv')
# set crs
crs = crs(conus_forestgroup)2.1.3 EPA level-3 Ecoregions
Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies
# level 3 ecoregions
l3eco <- st_read('data/ecoregion/us_eco_l3.shp') %>%
st_transform(., crs=crs)
# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>%
filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))2.1.4 Mapping
2.1.4.1 Ecoregions
# mapview
palette <- brewer.pal(18,"YlGn")
palette[1] <- rgb(255, 255, 255, maxColorValue=255, alpha=1)
mapview(eco_select,na.color=palette[1],zcol = "US_L3NAME",layer.name = "Level-3 Ecoregion", legend=TRUE)2.1.4.2 Forest Type Groups
# convert raster values to factors
forestgroup_eco <- crop(conus_forestgroup,eco_select) %>%
mask(.,eco_select)
forestgroup_eco[forestgroup_eco %in% c(120,180,240,300,320, 360, 370, 400, 500,700,900,920,950)] <- 1
# add a labels for forest type code
group_levels <- levels(forestgroup_eco %>% as.factor())[[1]]
group_levels[["Forest Type"]] <- c("Unforested","Other","200: Douglas-fir","220: Ponderosa Pine","260: Fir/Spruce/Mountain Hemlock","280: Lodgepole Pine")
# group_levels[["forest_type"]] <- c("Unforested","120: Spruce/Fir","180: Pinyon/Juniper","200: Douglas-fir","220: Ponderosa Pine","240: Western White Pine","260: Fir/Spruce/Mountain Hemlock","280: Lodgepole Pine","300: Hemlock/Sitka Spruce","320: Western Larch","360: Other Western Softwood","370: California Mixed Conifer","400: Oak/Pine","500: Oak/Hickory","700: Elm/Ash/Cottonwood","900: Aspen/Birch","920: Western Oak")
levels(forestgroup_eco) <- group_levels
# mapview
mapview(forestgroup_eco, col.regions=palette,na.color=palette[1],legend=TRUE,layer.name = "Forest Type")2.2 Define Fire Parameters
2.2.1 Monitoring Trends in Burn Severity (MTBS) Dataset
# import mtbs fire perimeters
mtbs_full <- st_read('data/mtbs/mtbs_perims_DD.shp') %>%
st_transform(., crs=crs)
# filter mtbs data to area and timepoints of interest
mtbs_select <- mtbs_full %>%
mutate(state = str_sub(Event_ID,0,2),
year = year(as.Date(Ig_Date))) %>%
filter(state %in% c("WA","ID","MT","WY","SD"),
between(Ig_Date, as.Date('1988-01-1'), as.Date('1991-12-31'))) 2.2.2 Group Adjacent Fires
# function to group adjoining fire polygons to ensure contiguous high-severity patches across MTBS events
group_fires <- function(mtbs_year) {
# join the polygons with themselves, and remove those that do not join with any besides themselves
combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>%
drop_na(Event_ID.y)%>%
dplyr::select(Event_ID.x,Event_ID.y)
if(nrow(combined)>=1){ # if there are overlaps for this years fires...
# partition data into that that has overlap, and that that does not
overlap <- mtbs_year %>%
filter(Event_ID %in% combined$Event_ID.x)
no_overlap <- mtbs_year %>%
filter(!(Event_ID %in% combined$Event_ID.x))
print(paste0("there are ",nrow(overlap)," overlapping polygons"))
# join all overlapping features, and buffer to ensure proper grouping
overlap_union <- st_union(overlap) %>%
st_buffer(190)
# break apart the joined polygons into their individual groups
groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
mutate(year = mean(mtbs_year$year),
Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
rename(geometry = x)
print(paste0("polygons formed into ",nrow(groups)," groups"))
# join back with original dataset to return to unbuffered geometry
grouped_overlap <- st_join(overlap,groups,left=TRUE)
# arrange by the new grouping
joined_overlap_groups <- grouped_overlap %>%
group_by(Fire_ID) %>%
tally()%>%
st_buffer(1) %>%
dplyr::select(Fire_ID) %>%
mutate(year = mean(mtbs_year$year))
# add new ID to the freestanding polygons
no_overlap_groups <- no_overlap %>%
mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
dplyr::select(Fire_ID,year)
# join the new grouped overlap and the polygons without overlap
fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
return(fires_export)
} else { # if there are no overlaps for this year...
print("no overlapping polygons")
fires_export <- mtbs_year %>%
mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
dplyr::select(Fire_ID,year)
return(fires_export)
}
}# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>% filter(year == 1988))## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"
fires_89 <- group_fires(mtbs_select %>% filter(year == 1989))## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_90 <- group_fires(mtbs_select %>% filter(year == 1990))## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_91 <- group_fires(mtbs_select %>% filter(year == 1991))## [1] "no overlapping polygons"
# join each fire year
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
mutate(area_ha = as.numeric(st_area(geometry))/10000,
area_acres = area_ha*2.471)2.2.3 Select Fires by Ecoregion and Forest Type
# assign ecoregion and proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>%
left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08,
fun = function(value, coverage_fraction) {
data.frame(value = value,
frac = coverage_fraction / sum(coverage_fraction)) %>%
group_by(value) %>%
summarize(freq = sum(frac), .groups = 'drop') %>%
pivot_wider(names_from = 'value',
names_prefix = 'freq_',
values_from = 'freq')}) %>%
mutate(across(starts_with('freq'), replace_na, 0)))
# remove unnecessary columns, cleanup names
# filter to ensure fire polygons are at least 25% type of interest
fires <- fires_join %>%
dplyr::select("Fire_ID","year","area_ha","area_acres","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>%
rename("ecoregion" = "US_L3NAME",
"freq_df"="freq_200",
"freq_pp"="freq_220",
"freq_fs"="freq_260",
"freq_lpp"="freq_280") %>%
mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
freq_forested = 1- freq_0,
freq_ideal = freq_df+freq_fs+freq_lpp)%>%
mutate(across(starts_with('freq'), round,2))%>%
filter(freq_ideal > 0.25)2.2.4 Select Fires by Burn Severity
# import all mtbs rasters via a list
rastlist <- list.files(path = "data/mtbs", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,22,25))
# create empty dataframe
severity_list <- list()
# loop through mtbs mosasics for 1988-1991
# extract mtbs burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
mtbs_year <- allrasters[[i]]
fire_year <- filter(fires, year==str_sub(i,2,5))
raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
names(raster_extract) <- fire_year$Fire_ID
output_select <- bind_rows(raster_extract, .id = "Fire_ID")%>%
group_by(Fire_ID , value) %>%
summarize(total_area = sum(coverage_area)) %>%
group_by(Fire_ID) %>%
mutate(proportion = total_area/sum(total_area))%>%
dplyr::select("Fire_ID","value","proportion") %>%
spread(.,key="value",value = "proportion")
severity_list[[i]] <- output_select
}
# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list)
# join burn severity % to fires polygons
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="Fire_ID")%>%
rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>%
dplyr::select(- "NaN",-"regrowth",-"error") %>%
mutate(highsev_acres = area_acres*highsev)%>%
filter(highsev_acres > 500)2.2.5 Clean Up Dataset
# get the most common forest type within each polygon
fires_select <- fires_severity %>%
left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08))
fires_select$mode <- as.factor(fires_select$mode)
fires_select <- fires_select %>%
mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
mode==220 ~ "Ponderosa",
mode==260 ~ "Fir-Spruce",
mode==280 ~ "Lodegepole Pine",
TRUE ~ "Other"),
Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year))# join the grouped fires back to original mtbs boundaries
fires_mtbs <- st_join(mtbs_select,fires_select,left=FALSE,largest=TRUE) %>%
filter(year.x==year.y)%>%
dplyr::select("Event_ID","Incid_Name","Fire_ID","Ig_Date","year.y","state","BurnBndAc","ecoregion") %>%
rename(year= year.y)2.3 Final Fire Dataset
2.3.1 Dataset Overview
full_dataset <- fires_select %>%
st_drop_geometry() %>%
dplyr::select(Fire_ID,year,ecoregion,fire_foresttype,area_acres,highsev) %>%
mutate(highsev = round(highsev,2),
area_acres = round(area_acres,0))
kable(full_dataset,
align = 'c',
padding = 1,
col.names = c("Fire ID", "Year", "Ecoregion", "Majority Forest Type","Area (acres)", "High Severity %"),
caption = "High-Severity Conifer-Dominated Fires 1988-1991")| Fire ID | Year | Ecoregion | Majority Forest Type | Area (acres) | High Severity % |
|---|---|---|---|---|---|
| Fire_1_1988 | 1988 | Middle Rockies | Fir-Spruce | 342005 | 0.27 |
| Fire_2_1988 | 1988 | Middle Rockies | Lodegepole Pine | 777690 | 0.28 |
| Fire_3_1988 | 1988 | Middle Rockies | Lodegepole Pine | 448911 | 0.21 |
| Fire_4_1988 | 1988 | Idaho Batholith | Douglas-Fir | 5651 | 0.16 |
| Fire_5_1988 | 1988 | Idaho Batholith | Fir-Spruce | 11945 | 0.23 |
| Fire_6_1988 | 1988 | Idaho Batholith | Douglas-Fir | 50666 | 0.07 |
| Fire_7_1988 | 1988 | Middle Rockies | Lodegepole Pine | 167870 | 0.50 |
| Fire_8_1988 | 1988 | Idaho Batholith | Douglas-Fir | 25889 | 0.02 |
| Fire_9_1988 | 1988 | Idaho Batholith | Douglas-Fir | 8312 | 0.17 |
| Fire_10_1988 | 1988 | Canadian Rockies | Lodegepole Pine | 42492 | 0.52 |
| Fire_11_1988 | 1988 | Idaho Batholith | Fir-Spruce | 45075 | 0.43 |
| Fire_12_1988 | 1988 | Middle Rockies | Lodegepole Pine | 5633 | 0.24 |
| Fire_13_1988 | 1988 | Idaho Batholith | Douglas-Fir | 4962 | 0.25 |
| Fire_14_1988 | 1988 | Middle Rockies | Other | 35864 | 0.46 |
| Fire_15_1988 | 1988 | Idaho Batholith | Fir-Spruce | 19499 | 0.29 |
| Fire_16_1988 | 1988 | Idaho Batholith | Fir-Spruce | 5626 | 0.09 |
| Fire_17_1988 | 1988 | Idaho Batholith | Fir-Spruce | 12746 | 0.40 |
| Fire_18_1988 | 1988 | Middle Rockies | Other | 13108 | 0.48 |
| Fire_19_1988 | 1988 | Idaho Batholith | Fir-Spruce | 7241 | 0.27 |
| Fire_20_1988 | 1988 | Idaho Batholith | Fir-Spruce | 6559 | 0.13 |
| Fire_21_1988 | 1988 | Middle Rockies | Fir-Spruce | 3113 | 0.17 |
| Fire_22_1988 | 1988 | Middle Rockies | Other | 29233 | 0.16 |
| Fire_23_1988 | 1988 | Middle Rockies | Other | 1363 | 0.41 |
| Fire_24_1988 | 1988 | Idaho Batholith | Douglas-Fir | 48136 | 0.31 |
| Fire_25_1988 | 1988 | Northern Rockies | Douglas-Fir | 8089 | 0.25 |
| Fire_26_1988 | 1988 | Middle Rockies | Douglas-Fir | 8588 | 0.56 |
| Fire_27_1988 | 1988 | Northern Rockies | Douglas-Fir | 11403 | 0.30 |
| Fire_28_1988 | 1988 | Northern Rockies | Douglas-Fir | 21854 | 0.25 |
| Fire_29_1988 | 1988 | Canadian Rockies | Douglas-Fir | 33844 | 0.27 |
| Fire_30_1988 | 1988 | Northern Rockies | Douglas-Fir | 1909 | 0.31 |
| Fire_31_1988 | 1988 | Middle Rockies | Other | 6282 | 0.47 |
| Fire_32_1989 | 1989 | Idaho Batholith | Fir-Spruce | 13334 | 0.15 |
| Fire_33_1989 | 1989 | Middle Rockies | Lodegepole Pine | 3298 | 0.37 |
| Fire_34_1989 | 1989 | Idaho Batholith | Douglas-Fir | 4928 | 0.38 |
| Fire_35_1989 | 1989 | Idaho Batholith | Douglas-Fir | 47680 | 0.02 |
| Fire_36_1989 | 1989 | Idaho Batholith | Fir-Spruce | 2486 | 0.30 |
| Fire_37_1989 | 1989 | Idaho Batholith | Fir-Spruce | 5566 | 0.30 |
| Fire_38_1989 | 1989 | Idaho Batholith | Fir-Spruce | 7443 | 0.28 |
| Fire_39_1989 | 1989 | Idaho Batholith | Fir-Spruce | 6786 | 0.26 |
| Fire_40_1989 | 1989 | Idaho Batholith | Douglas-Fir | 8733 | 0.12 |
| Fire_41_1989 | 1989 | Idaho Batholith | Fir-Spruce | 1615 | 0.49 |
| Fire_42_1989 | 1989 | Idaho Batholith | Lodegepole Pine | 2488 | 0.33 |
| Fire_43_1989 | 1989 | Idaho Batholith | Fir-Spruce | 3081 | 0.27 |
| Fire_44_1990 | 1990 | Middle Rockies | Douglas-Fir | 2535 | 0.45 |
| Fire_45_1990 | 1990 | Idaho Batholith | Fir-Spruce | 3418 | 0.53 |
| Fire_46_1990 | 1990 | Idaho Batholith | Douglas-Fir | 2249 | 0.49 |
| Fire_47_1990 | 1990 | Middle Rockies | Lodegepole Pine | 2763 | 0.41 |
| Fire_48_1990 | 1990 | Middle Rockies | Other | 13461 | 0.19 |
| Fire_49_1991 | 1991 | Middle Rockies | Douglas-Fir | 6978 | 0.31 |
| Fire_50_1991 | 1991 | Middle Rockies | Douglas-Fir | 3097 | 0.18 |
| Fire_51_1991 | 1991 | Middle Rockies | Fir-Spruce | 6995 | 0.14 |
| Fire_52_1991 | 1991 | Idaho Batholith | Douglas-Fir | 1186 | 0.55 |
| Fire_53_1991 | 1991 | Northern Rockies | Fir-Spruce | 7095 | 0.39 |
| Fire_54_1991 | 1991 | Northern Rockies | Fir-Spruce | 2478 | 0.51 |